/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package com.compomics.spectral.utilities.filtering.wavelet.daubechies.tresholding;

import org.apache.commons.math.stat.descriptive.DescriptiveStatistics;
import org.apache.commons.math.stat.descriptive.SummaryStatistics;

/**
 *
 * @author Kenneth
 */
public class ThresholdGenerator {
    /*The threshold
     (N being the signal length, σ being the noise variance) is well known in wavelet literature as the Universal threshold. 
     */

    public static double getThresHold(double[] intensities) {
        double[] tempIntensities = intensities.clone();
        double threshold = 0;
        double dev = 0;
        double avIntensity = 0;
        DescriptiveStatistics sumStats = new DescriptiveStatistics();

        //load summary
        for (double anIntent : tempIntensities) {
            sumStats.addValue(anIntent);
        }
        avIntensity = sumStats.getPercentile(50);
        dev = sumStats.getStandardDeviation();
        sumStats.clear();
        //remove all outliers?
        for (double anIntent : tempIntensities) {
            if (anIntent < (avIntensity + (3 * dev)) && anIntent > (avIntensity - (3 * dev))) {
                //Then this is NOT an outlier ---> only need noise
                sumStats.addValue(anIntent);
            }
        }
        //genuine noise deviation = 
        dev = sumStats.getStandardDeviation();
        int N = intensities.length;

        threshold = Math.abs(dev * Math.sqrt(2 * Math.log(N)));

        return threshold;
    }
}
